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We show that semiclassical formulas such as the Gutzwiller trace formula can be implemented 
on a quantum computer more efficiently than on a classical device. We give explicit quantum 
algorithms which yield quantum observables from classical trajectories, and which alternatively test 
the semiclassical approximation by computing classical actions from quantum evolution. The gain 
over classical computation is in general quadratic, and can be larger in some specific cases. 
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I. INTRODUCTION 



It is now widely recognized that the principles of quan- 
tum mechanics allow to realize new computational de- 
vices which can be more efficient than their classical 
counterparts [J |3) S 3 • Quantum algorithms have been 
proposed which take advantage of the quantum mechan- 
ical properties of these devices to perform specific tasks 
faster than on a classical computer. The most famous 
such algorithm is due to Shor [H and factors large in- 
tegers exponentially faster than any known classical al- 
gorithm. Another algorithm, for which the gain is onl 
quadratic, enables to search an unsorted database 
Efforts have been devoted also to using such quantum 
computers to simulate the behavior of complex physical 
systems, a task of much practical interest. Algorithms 
have been set up enabling to simulate certain quantum 
mechanical systems efficiently 0, H, H, [13, Eli: as was 
originally envisioned by Feynman. However, as quantum 
algorithms use procedures different from classical algo- 
rithms, it is by no means obvious which problems can be 
sped up by using a quantum computer. It is therefore 
important to precisely specify the class of problems that 
can be solved efficiently on a quantum computer, espe- 
cially among problems which have been implemented by 
scientists on classical devices because of their practical 
interest. 

On classical computers, a great deal of activity in the 
past decades has been devoted to the numerical imple- 
mentation of semiclassical formulas. Such formulas ap- 
proximate quantum mechanics through classical quanti- 
ties, and have been used since the beginning of quantum 
mechanics. Although they have been much studied, their 
application to practical computation of quantum observ- 
ables is often hampered by the exponential proliferation 
of classical orbits involved when the system is chaotic. 
Semiclassical formulas enable to approximate the exact 
quantum mechanics for small h, and give an insight into 
the relationship between classical and quantum mechan- 
ics. For integrable systems with n degrees of freedom, 
classical dynamics takes place on rt-dimensional tori in 
the 2n-dimensional phase space. In this case, semiclassi- 
cal formulas quantize individual tori. They are relatively 
straightforward to implement and have been constructed 



and used early in the development of quantum mechan- 
ics. In contrast, for chaotic systems this quantization of 
tori is not valid, as pointed out by Einstein as early as in 
1917 [12], and individual wavefunctions cannot be built 
from a single classical structure. As a substitute, various 
formulas have been constructed, which express the quan- 
tum quantities in terms of an (infinite) set of classical 
orbits. The most famous such formula is the Gutzwiller 
trace formula [l3| . where the quantum density of states 
d-{E) = X^n ^{E~En) (where En are the energy levels) is 
written as a function of all classical periodic orbits of the 
system. It has the general form d{E) = Yip ^p^'''''^^^ i 
where the sum runs over all periodic orbits, Ap is re- 
lated to the stability of the orbit and ipp to its action. 
It can be viewed as a Fourier-type duality between the 
set of all eigenenergies of the system on the one hand 
and the set of all actions of periodic orbits on the other 
hand. Other formulas of the same kind give the quan- 
tum propagator G{x,x') in terms of all classical orbits 
from X to x' (Van Vleck formula) or scattering ampli- 
tudes in term of scattering orbits (Miller's formula) [l^ . 
Many works have implemented numerically such formu- 
las by truncating the sum over classical orbits (see e.g. 
[H [ll m, m, [la, mi), e.g. to obtain the semiclassical 
spectrum, but because of the exponential proliferation of 
classical orbits typical of chaotic systems only a few semi- 
classical eigenvalues can be extracted. Several methods 
have been devised to reduce the number of orbits enter- 
ing the sum [2^ . [23l . [23 | , but they all require summing up 
contributions from a still exponential number of orbits. 

In this paper, we study the implementation of semi- 
classical formulas on quantum computers. We show that 
for certain dynamical systems, such formulas can be com- 
puted more efficiently on a quantum computer than on a 
classical device. From the quantum information point of 
view, this gives new examples of algorithms where a gain 
can be reached compared to classical algorithms. From 
the point of view of quantum chaos, this would enable 
these formula to become more practical on a quantum 
computer if such a device becomes available, and thus to 
explore the quantum-classical correspondence in regimes 
which are difficult to reach on a classical computer. The 
paper is organized as follows. In section II, we present 
in detail the most famous semiclassical formula which 
relates the density of states to classical periodic orbits 
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(Gutzwiller trace formula), in the specific case of quan- 
tum maps. We then discuss in section III a quantum 
algorithm which implements this semiclassical formula 
in the form where it is most difficult classically, i.e. sum- 
ming up classical orbits and extracting quantum observ- 
ables. In section IV we implement the same formula but 
in the reverse direction, i.e. using quantum observables 
to extract classical quantities. Our method can be con- 
sidered as a new way of extracting information from the 
quantum simulation of quantum systems. Indeed, while 
many quantum systems can be simulated efficiently on 
a quantum computer, a crucial point to get a complete 
algorithm and make the gain effective is to devise a read- 
out method once the simulation is performed. It has been 
shown that the gain over classical computation can de- 
pend critically on the observable measured at the end of 
the simulation [2^, [2^, [23| ■ In the present paper we show 
that in general we can expect a quadratic gain over classi- 
cal computation using the algorithms of sections III-IV, 
and that this gain can be quartic for some quantities. 
The original hope of this study was to use the quantum 
Fourier transform which is exponentially faster than the 
classical Fourier transform to ensure an exponential gain 
for this type of problem. It turned out that for most sys- 
tems counting the total number of gates involved shows 
that only a polynomial gain can be reached. However, in 
section V, we give an example of a related problem where 
exponential gain can be reached. 

II. SEMICLASSICAL TRACE FORMULA FOR 
QUANTUM MAPS 

Classical and quantum maps represent a particularly 
simple class of dynamical systems. Indeed, such systems, 
where one iteration of the map corresponds to a discrete 
time step, are easier to handle and yield simpler formulas. 
In what follows, we will restrict ourselves to such systems. 
This does not entail a major loss of generality, since it is 
known that Hamiltonian systems can in general be rep- 
resented by maps through the construction of Poincare 
surfaces of section [28j. Furthermore, most of the phe- 
nomena observed in more complicated systems can be re- 
produced in well-known models of quantum maps. This 
explains why many works on semiclassical formulas have 
used classical and quantum maps as testbeds. 

Here we consider two-dimensional maps on a toroidal 
phase space. Let us first give examples of well-known 
classical maps that we will use later on. A much studied 
instance is the family of cat maps [H, [2^, [s^, [3l|j i-G- 
linear automorphisms of the torus characterized by 2 x 2 

matrices of SL{2, Z). For a matrix M = ^ ^" ^ , the 

corresponding map is 

p ~ tiip + ti2q (modi) 

q = t2ip + t22q (modi), (1) 

where {p, q) are phase-space variables and bars denote 



new variables after one iteration of the map. 

Another well-known example is the baker's map f28j: 

{q,p) = (2g,|) forO<(7 < i 

{q,p) = (2g-l,^) fori<9<l. (2) 

Maps (II])-© are instances of strongly chaotic systems, 
with homogeneous exponential divergence of trajectories, 
positive Kolmogorov-Sinai entropy, and exponential pro- 
liferation of periodic orbits with the length. 

More generally, many classical maps can be written in 
the form 

p = p — kV'{q) 

q = q + Tp, (3) 

where the potential V{q) is a function of position. Such 
maps correspond to the integration over one period of 
a free rotator periodically kicked by a potential V{q). 
They include the standard map (the classical version of 
the kicked rotator) ^] for V{q) = cosq, or the saw- 
tooth map [l3| for V{q) — ~{p — vr)^/2. These maps 
display a wide range of different behaviors depending on 
the parameters. In particular, for the standard map the 
dynamics changes from close to integrability for small 
values of the parameter kT to fully developed chaos for 
large values of kT. 

The quantum version of the classical maps acts on a 
Hilbert space of dimension N corresponding to the in- 
verse of Planck's constant 2TTh. It is represented by an 
N X N matrix U I3p|. In the case of a cat map ([T|) , the 
quantization yields jsd [Blj 

where S{qi,q2) = {hiql - 2qiq2 + t22ql)/{2ti2) and the 
average is taken over all integers m. 

The quantized baker's map (ssj is even simpler. The 
evolution operator on a iV-dimensional space is given by 

where F„ is the N x N matrix with {Fn)kj = «^ 
(discrete Fourier transform). 

At last, maps of the form ^ yield, upon quantization, 
quantum maps of the form: 

U = f,-lTp'/2h^-lkV(q)/h_ ^g-j 

These evolution operators can be implemented effi- 
ciently on a quantum computer. This was shown for ([5]) 
in [3J] using the quantum Fourier transform instead of 
the classical one, and in [ol, fioj for maps of the form ([6]). 

One of the advantages of maps over generic systems is 
that some of the steps leading to the trace formula link- 
ing the spectrum to periodic orbits can be made exact. 
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Indeed, the spectral density for an N x N quantum map 
U with eigenphases 6*^., 1 < A; < A^, is given by 
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^ {e-^'hvU' + e^'hiU-') . (7) 



2n 2n ^ 

t=i 

This expression, obtained by Poisson summation for- 
mula, is exact and only depends on the traces of iterates 
of the quantum map. Similarly, one can express the co- 
efficients of the characteristic polynomial det (/ — xU) — 
^Ylik l^kx'' only in terms of traces of powers of U by using 
the recurrence relation 



1 ^' 

= -tE'^'=-**^"^*' /^"^l- 



(8) 



This relation can be easily proved by expanding dct(/ + 
zU) = exp tr log(/ -I- zU) into powers of z. Unitarity of 
the operator U implies the symmetry relation 



(3N-k - det{-U)pk- 



(9) 



Thanks to this resurgence relation the computation of 
trJ7* for t < N/2 suffices to calculate the characteristic 
polynomial. 

The semiclassical approximation of the spectrum can 
be obtained by calculating the coefficients of the char- 
acteristic polynomial ([8]) using semiclassical expressions 
for the traces. For large N the main contribution to 
trJ7* comes from periodic orbits. For a classical map 
(j) mapping the phase-space onto itself, a periodic orbit 
of length t is a fixed point of (/)*. It is given by a se- 
quence (poi 9o,Pi, 917 •■•jPtj 9t) of phase-space points such 
that {pi,qi) — 4>{pi-i,qi^i) for all z, 1 < i < t, and 
{Pt^Qt) — (Po^Qo)- If tp is the smallest integer such that 
{Ptp,qtp) = {po,qo), then tp divides t and the periodic or- 
bit is the repetition of r = t/tp times a primitive periodic 
orbit. A given primitive periodic orbit is characterized 
by its monodromy matrix Mp (which is the linearized 
version of the map (j) in the vicinity of the periodic or- 
bit), its action Sp = J2}'=i S{q^-i,qi) where S{q, q') is y g^ch that Ty = 1 if the image of Ri has a non-empty 



properties make difficult the practical use of semiclassi- 
cal formulas, which need enormous numbers of orbits to 
be accurate. As we will show in the next section, this 
task can be made easier on a quantum computer. The 
fact already mentioned that the quantum evolution op- 
erator of these maps can often be implemented efficiently 
on a quantum computer opens the way to the use of semi- 
classical formulas in the reverse direction, using quantum 
observables to infer results on classical quantities. This 
will be the subject of section IV. 



III. SPECTRUM FROM CLASSICAL 
QUANTITIES 

We first discuss an algorithm allowing to obtain semi- 
classically the set of eigenvalues of the quantum map, 
or equivalently the coefficients ^ of the characteristic 
polynomial of the map. 

In order to calculate the traces using (fTO|) we need to be 
able to characterize periodic orbits of the classical map. 
There are instances of systems where this task is very 
easy. For instance for cat maps H]), the iterates of the 
classical map can be calculated analytically, and there- 
fore periodic orbits are entirely characterized. This is also 
the case for perturbed cat maps which are Anosov maps 
of the form (j) — (f>Q o Xe, where 0o is a cat map and Yp 
is a perturbation close to the identity. It was shown [36[ 
that for sufficiently weak perturbations orbits of Anosov 
maps remain topologically conjugate to periodic orbits of 
the unperturbed cat map . Thus periodic orbits can be 
described completely [37|. More generically we will con- 
sider systems in which periodic orbits can be described 
by a symbolic dynamics associated with a finite Markov 
partition. That is, phase space can be partitioned into 
sets 1 < fc < TO, and intersections of the images of 
the Rk under the (forward and backward) iterates of the 
classical map define finer and finer partitions so that at 
infinity the intersections contain either no point or a sin- 
gle one [ll]. Thus a given (infinite) sequence of labels 
corresponds to at most one point of phase space. The 
mapping rules between the Rk under one iteration of the 
map can be summarized in a to x to transition matrix 



the classical action from q to q' , and its Maslov index Vp 
The semiclassical expansion of trC/* reads 



f gir(Sp/h-VpTr/2) 

''^'""*^,§jdet(J-M;)|V2' 



(10) 



where the sum runs over the set Vt of all periodic orbits 
of length (number of time steps) t = rtp. The action, 
Maslov index and monodromy matrix correspond to the 
associated primitive periodic orbit [s^ . 

For maps (IJl-Q, the classical dynamics displays some 
form of chaos, up to the strongest types with exponential 
divergence of nearby trajectories and exponential prolif- 
eration of periodic orbits with increasing length. Such 



intersection with Rj, and otherwise. This transition 
matrix sums up the grammar rules that discriminate be- 
tween allowed words and forbidden ones. There is a one- 
to-one correspondence between phase-space points and 
allowed symbolic sequences, and periodic orbits corre- 
spond to periodic sequences of symbols. 

Simple examples of quantum maps with symbolic dy- 
namics are perturbed cat maps, or the baker's map [39j . 
In the latter example, symbolic dynamics is described by 
only two symbols and 1, and all sequences of symbols 
are allowed. Dynamical systems such as the 3-disk ji^ 
or motion on surfaces with constant negative curvature 
[4l| also provide examples where a symbolic dynamics 
exists with all sequences allowed. In such examples. 
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periodic trajectories are in one-to-one correspondence 
with periodic strings of and 1. 



(as in Fig. [T]) but also by qubits of register B. Steps I 
and II are polynomial in t^^ax- 



We will now sketch the steps of a quantum algo- 
rithm allowing to compute the semiclassical traces 
()10|1 in a parallel way. To simplify notation, for 
each periodic orbit p of length t = rtp we define the 
amplitude Ap — ip/|det(J — Mp)|^/^ and the phase 
(j)p — r{Sp/h~ VpT^I^). Thus we have to calculate the 
quantities Tf = X]p -^p^*'^''- consider a system 

whose symbolic dynamics is described by a finite Markov 
partition. For simplicity we assume that the partition 
consists of only two sets. Then only two symbols and 
1 are required (if there are more than two sets in the 
partition we code labels by binary strings). The trace 
formula \\S)\ will be truncated at imax, which means that 
only periodic orbits of length t < iniax will be considered. 
We distinguish five registers in the computational state. 
Register A will hold the lengths t, < t < t^^x of the 
periodic orbits. It requires riA = log2(imax) qubits. Reg- 
ister B will hold, on its last t qubits, the 2* codewords 
corresponding to a given orbit length t. This register 
has to contain tib = imax qubits. Register C is used 
to store the phases (j)p, and register D is used for the 
"tuning" of the amplitudes Ap associated to each orbit. 
Additional registers will serve as workspace. We will 
make use of following one-qubit operations: rotation of 



and phase shifts 



the A;th qubit Rk{0) = exp 
Pk{9) — exp (^—iOai'^''^. The steps are as follows. 

Step I: Let A be such that the number of allowed 
codewords scales exponentially with t as exp(Ai). 
We define A such that for each t, the amplitude 
Ap = tp/\det{I — Mp)\^/^ of each periodic orbit p of 
length t = rtp is upper bounded by exp(— Ai). The 
amplitude Ap in many cases will be actually close 
to exp(— At). Let us set /i = A — A/2. We define 

angles Ok e [0,n/2] by cos6lfc = + e-^''^" . Ap- 

plying riA rotations Rk{dk) to register A of the initial 
state gives (up to a normalization factor) the state 

Step II: All allowed codewords are generated on reg- 
ister B. In the simplest case where there is no grammar 
rule one wants, for each value t on register A, to put 
the last t qubits of B into a uniform superposition. This 
is performed by applying, for each t, t Hadamard gates 
controlled by register A on the last t qubits of B (see 
Fig- HI) ■ This gives (up to normalization) the state 



EE' 

t p 



-At 



\t)A\p)B\0)c\0)D, 



(11) 



where the second sum runs over all codewords p, 
< p < 2* — 1. If there is a finite number of grammar 
rules the allowed codewords can be generated by replac- 
ing the Hadamard gates by rotations Rk- For each value 
of t these rotations are controlled not only by register A 



register A 



register i? 



H — H 



H 



H — 



H 



H 



FIG. 1: Circuit for step II and two-letter symbolic dynamics. 
Register A codes for lengths t, < t < 3 on two qubits. 
The Hadamard gates are controlled by the values of t, and on 
register B the state \t)\0) becomes 2"*/^ Eilo^ I*)!*)- 

Step III: From each codeword p of length t it 
is possible to recover the phase-space coordinates 
{po,qo,Pi,qi, ■■■,Pt,qt) of the trajectory coded by this 
codeword, as well as the characteristics of this trajec- 
tory: action Sp, monodromy matrix Mp, Maslov index 
Vp. These quantities can be calculated in a parallel 
way by classical operations implemented on the quan- 
tum workspace registers, as has been done classically in 
many systems [M [13, [11, [H, [H [Hj . The values of 0p 
and In Ap + At are then calculated and written on regis- 
ters C and D. For the kind of systems considered here, 
this step is polynomial in tmax as the operations are per- 
formed in parallel. After erasing intermediate steps we 
get a state 

J2 E e-'''\t)A\p)B\^p)c\ InAp + At)D. (12) 
t peVt 

Step IV: As in step I we use rotations Rk{Ok) to trans- 
fer the value stored in register D into an exponential 
prefactor exp(ln^p -I- M). The angles 9k are now given 

by cos^fe 



1/ V 1 -I- e~2'*2 , where the constant k sets 
the precision that one wants to achieve on the prefactor. 
For each value of t and p, register D is the sum of two 
orthogonal components \Q)d and |'0p)d- As the prefac- 
tor e^^* in (|12p is meant to yield a rough estimate of 
the amplitudes Ap, it can be expected that the quanti- 
ties In Ap + At are small and that the relative weight of 
(ipplipp) is small. 

By controlled phase shifts on register C the states are 
then multiplied by the phase factor e*"^', and step III is 
run backwards to erase register C. This yields 

J2 Ape^^^\t)A\p)B\0)cmD + \^Pp)d) (13) 
* peVt 

Step V: In order to get the semiclassical traces Tt — 
^pApC^'^'', we perform imax Quantum Fourier Trans- 
forms (QFT) on register B. Each QFT corresponds to a 
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given value of t and operates on the last t qubits of reg- 
ister B. That is, the gates of the QFT are controlled by 
register A (as in step II, see Fig. [T|). This yields a state 



2'-l 

E E E Ae'^'-e-''^'''^''\t)A\k)B\0)cmD + \^p)D). 

t k=0 pGVt 

(14) 

Step VI: In the amphtude of the |A; = 0)_b term 
corresponds to the semiclassical traces Tt- Therefore we 
now just have to perform a quantum search of 1 0) b | 0)_d in 
()14|) . This is done by amplitude amplification performed 
on registers B and D. This process, which is the slowest 
part of our algorithm, requires 0(2*™^'=/^) operations (as 
{ipp\ipp) is small, the search on register D is expected to 
contribute only a prefactor). It brings the state into 
a state 



J2 E Ae'''^\t)A\0)B\0)c\0)D ^Y.^t\t)A\0)B\0)c\0)D. 
t p&pt t 

(15) 

Quantum state tomography then gives the relative values 
of all semiclassical traces Tt . The knowledge of ti (easily 
computed classically) allows to obtain the absolute 
values of the r*, and thus the characteristic polynomial. 
Because of the symmetry relation (|9]) only traces up to 
^max — N/2 are required. Therefore the cost of our quan- 
tum algorithm (which is essentially the cost of amplitude 
amplification in step VI) is 0{2^/^). This is to be com- 
pared with the classical cost of 0(2^/^) required for the 
calculation of the semiclassical characteristic polynomial. 

As already mentioned our algorithm aims at estimat- 
ing the accuracy of the semiclassical approximation. Ob- 
viously, the cost of calculating the exact characteristic 
polynomial, with a scaling in 0{N'^), is far less. Thus for 
systems where the trace formula is exact, the result of 
the semiclassical sum should only yield with much more 
efforts the same result as the exact diagonalization. For 
instance for cat maps the exact equality tr?7* = Tt holds 
in Eq. pTI)) . and thus cat maps are not suited to study- 
ing discrepancies between exact and semiclassical energy 
levels if the full semiclassical sum is used. There are 
however instances of systems for which characterization 
of classical periodic orbits remains easy while the traces 
obtained through the trace formula (|10p are truly approx- 
imations, such as e. g. the perturbed cat maps described 
above. In such cases, our algorithm yields the semiclassi- 
cal spectrum with quadratic efficiency compared to clas- 
sical computation. Besides, even when the trace formula 
is exact, its truncation is not, and therefore its implemen- 
tation has some interest and has been done classically in 
[H, [13, IHl ■ Indeed, it enables to understand the conver- 
gence properties of the sum over periodic orbits in ([7]). 



IV. CLASSICAL ORBITS FROM QUANTUM 
OPERATOR 

Another way of estimating the accuracy of the semi- 
classical approximation is to calculate how well the clas- 
sical actions of the periodic orbits are reproduced when 
calculated from the spectrum through the trace formula. 
In the semiclassical approximation the trace of the it- 
erates of the quantum operator can be written as a sum 
over periodic orbits. This sum can be put under the form 
Tt = Ae^"^^" ~ trC/* (see Eq. The actions 

Sp calculated from the quantum spectrum through the 
semiclassical formula (|10p are obtained by performing a 
Fast Fourier Transform (EFT) on the set of traces tr[/* 
calculated for all matrix sizes < < A^max- The num- 
ber of traces A^max to evaluate depends on the precision 
required for the actions. 

We now discuss a quantum algorithm allowing 
to calculate each trace tr[/*, for any matrix size 
< A^ < Amax- Let m be the smallest integer such 
that A^ < 2™, and M = 2™. We distinguish three 
registers in the state vector on which the computation is 
performed. Register A stores the lengths t of the orbits, 
<t < tmax, on riA = log^inax qubits; here tmax is some 
fixed integer specifying the highest period that we want 
to consider. The two other registers, each of length m, 
will store the computational basis vectors. Additional 
workspace registers will be used as well in the course of 
the computations. Starting from the state |0)yi|0)B|0)c, 
we perform the following steps. 

Step I: We first apply Hadamard gates on registers A, 
B and C to put them in an equal superposition of basis 
vectors. We obtain 

^ J2 \t)A\^)B\^)c. (16) 
t i=0 

What we want is in fact a sum running over a range 
0<i<A^— 1. To obtain this from (|16p we use an 
auxiliary qubit (register D) that is set to |0) if i — A' < 
and to |1) if i — A^ > 0. The relative weight of the state 

AT-l 

J2J2\*)^\')^\')o\o)d (17) 

t i=0 

is greater than 1/2. 

Step II: The N x N matrix [/* has to be applied to 
register B of each state |i)A|*)s|i)c- As an illustration 
we focus on operators of the type ([6]). It was shown in 
Pi that for such maps one iteration can be implemented 
efficiently for a fixed matrix of size a power of 2. The 
algorithm consists in using QFTs to shift back and forth 
between p and q representation, while the operators 
g'/(p) and e*^^'''' are applied in the basis where they 
are diagonal by multiplication of basis vectors by a 
phase. For TV ^ 2™ the simulation of the quantum map 
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involves a QFT on vectors of size not a power of 2. Such 
a procedure was proposed in [43| for any fixed vector 
size iV. The simulation of U can therefore be done 
efficiently. The simulation of C/* can be done sequen- 
tially, controlled by the qubits of register A (as in Fig.[T]). 

Step III: The state ([T7| is now transformed into: 

N-l N~l 

J2 E 1*)^ \^)c\0)D -EE Ul\t)A\J)B\^)c\0) 

t i t i,j—0 

(18) 

By amplitude amplification on registers _B, C, D we select 
vectors with \j)B = \i)c and \0)d, leading to 



E E U!Jt)A\^)B\^)c\0)D■ 



(19) 



i=0 



After erasing register C we perform a QFT on register 
B. As in section IIIIl we use amplitude amplifica- 
tion to select the state \0)b, whose amplitude is 
Y^i^UVM = trU^/^/M. This is the slowest step 
in our computation. For chaotic systems the matrix 
elements C/*,- for N x N matrices are of order I/Vn and 
the traces tr?7* are expected to be of order 1. Thus each 
amplitude amplification has a cost 0{^/N) and step III 
requires N Grover iterations in total. For integrablc 
systems the traces are of order \/lV, and therefore only 
one of the amplitude amplifications is needed, requiring 
TV Grover iterations in total for step III. 



Step IV: We are now in the state 

EtrC/*|i)A|0)s|0)c|0)j 



(20) 



The relative values of the traces for different values 
of t are obtained by quantum state tomography. The 
traces themselves are then deduced from the classical 
calculation of trJ7, requiring 0{N) classical operations. 

The algorithm requires the calculation of N^e^x traces 
trf/*, with < < A^max- Thus the cost of the quantum 
algorithm is of order A^^^x- Classically, we need to com- 
pute all A^niax traces. Except for trf/jv this would need 
0{N'^) classical operations if the map is of the type ([6]), 
and up to 0{N^) in the general case where diagonaliza- 
tion of the operator is required. Thus the classical cost 
is of order A^^^^ to A^j^^^^ operations. Thus in both cases 
the quantum algorithm outperforms classical computa- 
tion, albeit polynomially. 

We note that if one is interested in distinguishing inte- 
grablc and chaotic systems via the form factor as in the 
algorithm proposed in fz^. then one needs only to be 
able to distinguish traces of order ^/N (integrable case) 
from ones of order 1 (chaotic case), for a specific value 
oi N. In this case one can stop at step III and check 
that ^/N Grover iterations are enough to get to the state 



|0), in which case one concludes that the system is inte- 
grable, or not enough, in which case one concludes that 
the system is chaotic. Our algorithm then only needs 
0{\/N) quantum operations instead of 0{N'^) classical 
operations, an improvement from the quadratic gain in 
[201 . One can also compute exactly the trace (stopping 
at step IV), and compute the form factor for small with 
a quadratic improvement compared to classical compu- 
tation. 

D- 

V. EXPONENTIAL SPEED-UP BY PHASE 
ESTIMATION 

The preceding processes can be applied to many phys- 
ical systems and yield a polynomial speed-up compared 
to classical computation. However there exist systems 
where a larger (up to exponential) gain might be ob- 
tained, following a different strategy based on phase es- 
timation. This method [1, |42| enables to obtain an 
eigenvalue of a given operator U by applying condition- 
ally iterates of U to an eigenvector \^)\ this gives the 
state which, by Fourier transform on the 

first register, gives where exp(z6') is the eigen- 

value corresponding to \^). If Y^) is not an eigenvec- 
tor but some randomly chosen state, the same process 
leads to "^j l^i) l^j) where are eigenvectors and 
\^) — X^j '^jl^j)- To be efficient, this method criti- 
cally requires not only that U can be efficiently imple- 
mented, but also that exponential iterates of U can be 
implemented with polynomial number of quantum gates, 
a much more stringent requirement. In the case of the 
quantum cat map, this method is efficient and remark- 
ably enough can lead to classical quantities with expo- 
nential efficiency. 

It is known [30| that the nth iterate of the quantized 
cat map ([4]) coincides with the quantization of the classi- 
cal nth iterate. In [43l |. it was shown that one can simu- 
late the classical cat map efficiently on a quantum com- 
puter, while in (4^ . it was further shown that one can 
compute the classical nth iterate for exponentially large 
n with polynomial number of gates. Thus if one starts 
from a random vector j^*), one can compute |i)J7'|5') 
in polynomial number of gates for exponential i's and 
N] a quantum Fourier transform followed by a quantum 
measurement leads to the value \9j) of one eigenvalue of 
the quantum cat map. It is known [sij that these eigen- 
values are very constrained, being of the form 



27rj + (j){N) 
^I(AO ' 



(21) 



where n{N) is the quantum period function, that is the 
smallest integer such that 



JjniN) ^ j^t^N) 



(22) 



and the phase (j>{N) can be calculated easily from the 
components of matrix L 31]. 
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Thus using this algorithm the quantum period func- 
tion can be obtained in polynomial time on a quan- 
tum computer. This quantity is related to the classical 
period function, which for each matrix L is the short- 
est integer g such that = I modA^. Indeed, the 
quantum period n{N) is also the smallest integer such 
that L"(^) = / (modTV) if N is odd, and such that 

= ( 1 « ^ if iV is even. The 



(mod2A^) 1 (modA^) J 
two functions in all cases differ by at most a factor of 
two [Slj , so knowing one of them enables to test and find 
easily the other one. The classical period function de- 
scribes the periodic orbits of the classical cat map. It 
has been shown in [31 that finding it is as complex as 
factorization of integers, and can nevertheless be realized 
on a quantum computer polynomially fast using a vari- 
ant of order-finding. The use of the quantum cat map 
enables to get this classical quantity by an equally effi- 
cient alternate quantum algorithm, showing that in this 
specific case classical quantities can be obtained through 
quantum mechanics with exponential efficiency compare 
to classical algorithms. 

VI. CONCLUSION 

In the studies above, we have shown that it is possible 
to implement semiclassical formulas on quantum com- 
puters, with a gain on efficiency over the implementation 
on a classical computer. The gain is in general polyno- 
mial, but in specific instances an exponential gain can be 
obtained for related problems. We mention again that 
the algorithms of Section IV can also be used to study 
quantum systems without reference to the semiclassical 
approximation, in the manner of [26| . with actually a 
larger gain. 

The algorithms of sections III and IV can be general- 
ized to a large class of systems. Indeed, to generalize sec- 
tion III one can use the tool of Poincare surface of section 
to transform systems with continuous time to discrete 



maps. For example, a popular system to study quantum 
chaos corresponds to biUiards, i.e. classically a parti- 
cle bouncing between walls, and quantum mechanically 
a wave function obeying Helmholtz equation with bound- 
ary conditions. In this case, a simple surface of section is 
represented by the boundary, the phase space coordinates 
being the curvilinear abscissa along the boundary and the 
angle that the outgoing trajectory makes with the vector 
normal to the boundary. An alternate possibility would 
be to stick with the continuous time dynamics and use the 
semiclassical formulas appropriate for this case. In both 
cases, it is important to be able to enumerate the clas- 
sical trajectories used in the semiclassical sums, which 
requires that a reasonably good symbolic dynamics can 
be constructed (e.g. with finite Markov partition). This 
is already the case for classical implementations of these 
formulas, which have all been performed in such cases. 
Additionally, the method exposed in section III is all the 
more efficient since the Lyapunov exponent of orbits is 
uniform. In the case where the stability of different orbits 
varies widely in different phase space regions, the quan- 
tum algorithm will become less efficient. Thus although 
strongly chaotic systems are the most difficult to treat by 
semiclassical formulas, they are probably the ones where 
the algorithms above will be the most efficient compared 
to classical algorithms. 

To generalize Section IV to systems with continuous 
time is probably possible, but would necessitate to first 
build an algorithm to simulate such systems on quantum 
computers. We think that once this is done, the main 
ideas of our algorithm in section IV should then be ap- 
plicable. 

The quantum algorithms presented here can be applied 
to a wide variety of systems. They show that in a domain 
where extensive numerical simulations have been used in 
the past decades, a quantum computer could significantly 
improve the speed of the calculations. 

We thank the French ANR (project INFOSYSQQ) and 
the IST-FET program of the EC (project EUROSQIP) 
for funding. 
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